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Computing tensor eigenvalues via homotopy methods 


Liping Chen* Lixing Han^* Liangmin Zhou^ 


Abstract 

We introduce the concept of mode-A: generalized eigenvalues and eigenvectors of 
a tensor and prove some properties of such eigenpairs. In particular, we derive an 
upper bound for the number of equivalence classes of generalized tensor eigenpairs 
using mixed volume. Based on this bound and the structures of tensor eigenvalue 
problems, we propose two homotopy continuation type algorithms to solve tensor 
eigenproblems. With proper implementation, these methods can find all equivalence 
classes of isolated generalized eigenpairs and some generalized eigenpairs contained 
in the positive dimensional components (if there are any). We also introduce an 
algorithm that combines a heuristic approach and a Newton homotopy method to 
extract real generalized eigenpairs from the found complex generalized eigenpairs. A 
MATLAB software package TenEig has been developed to implement these methods. 
Numerical results are presented to illustrate the effectiveness and efficiency of TenEig 
for computing complex or real generalized eigenpairs. 


Key words, tensors, mode-Zc eigenvalues, polynomial systems, homotopy continuation, 
TenEig. 
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1 Introduction 

Eigenvalues of tensors were first introduced by Lim [29] and Qi m in 2005. Since then, 
tensor eigenvalues have found applications in automatic control, statistical data analy¬ 
sis, diffusion tensor imaging, image authenticity verification, spectral hypergraph theory, 
and quantum entanglement, etc., see for example, [SKIIlEIlEZlESlIinillllllS] and the 
references therein. The tensor eigenvalue problem has become an important subject of 
numerical multilinear algebra. 

Various definitions of eigenvalues for tensors have been proposed in the literature, in¬ 
cluding E-eigenvalues and eigenvalues in the complex field, and Z-eigenvalues, H-eigenvalues, 
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and D-eigenvalues in the real field [2911371 HO] . In |7], Chang, Pearson, and Zhang intro¬ 
duced a notion of generalized eigenvalues for tensors that unifies several types of eigenval¬ 
ues. Recently this definition has been further generalized by Cui, Dai, and Nie m- 

Unlike the matrix eigenvalue problem, computing eigenvalues of the third or higher 
order tensors is a difficult problem m- Nonetheless, several algorithms which aim at com¬ 
puting one or some eigenvalues of a tensor have been developed recently. These algorithms 
are designed for tensors of certain type, such as entry-wise nonnegative or symmetric ten¬ 
sors. 

For nonnegative tensors, Ng, Qi, and Zhou m proposed a power-type method for 
computing the largest H-eigenvalue of a nonnegative tensor. Modified versions of the 
Ni-Qi-Zhou method have been proposed in [aiiisoiEi]. 

For real symmetric tensors, Hu, Huang, and Qi [20] proposed a sequential semidefinite 
programming method for computing extreme Z-eigenvalues. Hao, Cui, and Dai m pro¬ 
posed a sequential subspace projection method for a similar purpose. Kolda and Mayo 
[23] proposed a shifted power method (SSHOPM) for computing a Z-eigenvalue. They 
have improved SSHOPM in [2l| by updating the shift parameter adaptively. The result¬ 
ing method can be used to compute a real generalized eigenvalue. Han [15] proposed an 
unconstrained optimization method for computing a real generalized eigenvalue for even 
order real symmetric tensors. The methods in can find more eigenvalues of 

a symmetric tensor if they are run multiple times using different starting points. Re¬ 
cently, Cui, Dai, and Nie |T2| proposed a novel method for computing all real generalized 
eigenvalues. 

In this paper, we are concerned with computing all eigenpairs of a general real or 
complex tensor. As can be seen from the next section, finding eigenpairs of a tensor 
amounts to solving a system of polynomials. Naturally one would consider to use algebraic 
geometry methods such as the Grobner basis method and the resultant method [10] for 
this purpose. These methods can obtain symbolic solutions of a polynomial system, which 
are accurate. However, they are expensive in terms of computational cost and space. 
Moreover, they are difficult to parallelize. A class of numerical methods, the homotopy 
continuation methods, can overcome these shortcomings of the Grobner basis and the 
resultant methods. During the past few decades, significant advances have been made on 
homotopy continuation methods for polynomial systems, see for example, [31 12^1^I321IT3] . 
Recently, the homotopy techniques have been used to study tensor decomposition and 
perfect identification problems ('[T7]b 

In this paper we investigate computing complex eigenpairs of general tensors using 
homotopy continuation methods. One attractive feature of the homotopy continuation 
methods is that they can find all isolated solutions of polynomial systems and some so¬ 
lutions in the positive dimensional solution components. We propose two homotopy-type 
algorithms for computing complex eigenpairs of a tensor. These algorithms allow us to 
find all equivalence classes of isolated eigenpairs of a general tensor and some eigenpairs 
in positive dimensional eigenspaces (if there are any). We also present an algortihm com¬ 
bining a heuristic approach and a Newton homotopy method to compute real eigenpairs 
based on the found complex eigenpairs. Numerical examples show that our methods are 
effective and efficient. 
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This paper is organized as follows. In Section 2, we define mode-/c generalized eigen¬ 
values and eigenvectors which extend the matrix right eigenpairs and left eigenpairs to 
higher order tensors. Some properties of such eigenpairs are proved. An upper bound 
for the number of equivalence classes of generalized tensor eigenpairs using mixed volume 
is derived. In Section 3, we consider computing mode-A: generalized complex eigenpairs 
and present two algorithms. In Section 4, we introduce a method to compute real mode-fc 
generalized eigenpairs. Finally in Section 5, some numerical results are provided. 


2 Tensor eigenvalues and eigenvectors 

Let F = C or M be the complex field or the real field. Let m > 2, m' > 2, and n be positive 
integers. Denote the set of all mth-order, n-dimensional tensors on the field F by fI™’"'!. 
A tensor in is indexed as 

A = (Aqij-"*™)) 

where G F, for 1 < ' ,im < n. 

For X G C”, the tensor A defines a multilinear form 


Ax'^ = Y, 


A. 




5 ■ ■ ■ — 1 


( 2 . 1 ) 


For 1 < k < m, ^ is an n-vector whose jth entry is defined as 




ill'" tik — liik-\-l')"‘ lim — 1 


' ' ' ^ik-i^ik+l ' ' ' • 


( 2 . 2 ) 


When k = 1, the vector is denoted by 

A real tensor A G is positive definite if the multilinear form Ax"^ is positive 

for all x G M"'\{0}. A tensor A G fI^’A is symmetric if its entries are invariant 

under any permutations of their indices ii, ^2, •'' Am- 

We now introduce the following mode-/c generalized eigenvalue definition for a general 
tensor A. 


DEFINITION 2.1 Let A G F[™’A and B G Assume that Bx^' is not identically 

zero as a function of x. For 1 < k < m, if there exist a scalar A G C and a vector 
x G C"'\{0} such that 

• when m = m', 

AA)x^-^ = XBx^-\ (2.3) 

• when m A 

^ XBx'^'-^, Bx^' = 1, (2.4) 
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then A is called a mode-k B-eigenvalue of A and x a mode-k B-eigenvector associated with 
A. (A,x) is called a mode-k B-eigenpair of A. 

// A € M, X E M", then A is called a mode-k Br- eigenvalue of A and x a mode-k 
BR-eigenvector associated with X, and (A, x) a mode-k Bn-eigenpair of A. 

Denote the set of all mode-k B eigenvalues of A by asiA^^^). 

REMARK 2.1 Let (A, x) be a mode-/c 0-eigenpair of A. By (I2.3p or (I2.4p . (A,x) is 
a solution to = XBx^'~^. So is (A',x') with A' = and x' = tx for 

t E C\{0}. From this point of view, the solution space of = XBx'^ consists 

of different equivalence classes. We denote such an equivalence class by 

[(A, x)] := {(A', x') I A' = t"'A, x' = tx, t E C\{0}}. 

When m / m', taking arbitrary (A', x') E [(A, x)] and substituting x' = tx into Bx™ = 
1 in (12.4p yields t™ =1, which gives m' different values for t. This implies that the 
normalization Bx'^ = 1 in (j2.4p restricts us to choose m! representative solutions from 
each equivalence class. 

In our later discussions, we often choose only one representative from each equivalence 
class, and we often count the number of equivalence classes of mode-A: B-eigenpairs. 

REMARK 2.2 If only one representative is desirable from each equivalence class of 
eigenpairs, we can solve = ABx™' augmented with an additional linear equa¬ 

tion 

aiXi -I- 02X2 H-h OnXn -|- 6 = 0, (2.5) 

where ai,... ,an,b are random complex numbers. Then normalize the resulting solution 
to satisfy Bx"* = 1 in the case m A 

In the matrix case when m = m' = 2 and B = In (the nxn identity matrix), the mode-1 
eigenvectors are right eigenvectors and the mode-2 eigenvectors are left eigenvectors of A, 
and the mode-1 and mode-2 eigenvalues are the eigenvalues of matrix A, i.e., (Tb(A^^^) = 
<^b{A^^'^). However, when m > 3, = aB{A^^^) is generally not true when k A U 

unless A has a certain type of symmetry. The following example illustrates this situation. 

EXAMPLE 2.1 Consider the tensor A E whose entries are 

^111 = 1, Ai2i = 2, A 211 = 3, A 221 = 4, 

^112 = 5, A 122 = 6, A 212 = 7, A 222 = 0. 

Choose m' = 2 and B = I 2 (the 2x2 identity matrix). Note that in this case, if (A,x) is 
an B-eigenpair of A, so is (—A,—x). We follow [5], regarding (A,x) and (—A,—x) as the 
same eigenpair. Then 

^^(A^^)) = {0.4105,4.3820,9.8995}, 
asiA^'^'’) = {0.2851,4.3536,9.5652}, 
asiA^^^) = {0.2936,4.3007,9.4025}. 

Clearly, ct/j(A*^^^) A (Tb(A^^^) when k A I- 
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PROPOSITION 2.1 Suppose that (A, x) is a mode-k B-eigenpair, (/x, x) is a mode-l 
B-eigenpair of A, and when m = m!, Bx^ ^ 0. Then A = /x. 


Proof: Note that 


\Bx^ 

This immediately implies that A = /x, 


= Ax"^ = 
since Bx'^' 


PlBx'^'. 

= 1 when m 7 ^ m' and Bx'^ A 0 when 


vn = m!. 


□ 


Let A € For 1 < k < I < m, tensor G € is said to be the {k, 1) transpose 

of A if 

for all 1 < ii, • • • , Xm < rn. Denote the {k,l) transpose of A by A^^’^K We say that tensor 
A is {k, 1) partially symmetric if 

= A. 

PROPOSITION 2.2 Let A € and B € Assume that Bx"^' is not identi¬ 

cally zero as a function of x. Let k,l be integers such that 1 < k < I < m. Then 

• (A,x) is a mode-k B-eigenpair of A if and only if it is a mode-l B-eigenpair ofA^^’^^. 

• The sets of mode-k B-eigenpairs and mode-l B-eigenpairs are the same if A is {k,l) 
partially symmetric. 

The eigenvalues/eigenvectors defined in [3[T2l|33|40] are mode-l eigenvalues/eigenvectors. 
The tensors considered in these papers are primarily real symmetric tensors. For symmet¬ 
ric tensors, the sets of mode-A: ;B-eigenpairs and mode-l ;B-eigenpairs are the same for any 
k. Therefore, mode-l eigenvalues serve the purpose of those papers. On the other hand, 
nonsymmetric tensors arise from applications and theoretical studies, see, for example, 

[a El [la EH sa 08]. In [29], Lim defined mode-A: eigenvalues/eigenvectors for nonsym¬ 
metric real tensors A when B is the m'ih order identity tensor for some m' >2. Definition 
12. II considers more general A and B. 

As in mm, Definition 12.11 adapts a unified approach to define tensor eigenvalues. It 
covers various types of tensor eigenvalues introduced in the literature, including 

• If A, € m' = 2, and B is the identity matrix € R"’^"', the mode-I B- 

eigenpairs are the E-eigenpairs and the mode-I ;Bij-eigenpairs are the Z-eigenpairs 
defined in m, which satisfy 

Ax^~^ = Ax, x'^x = 1. (2.6) 

• If A G R[™'’"’1, m' = 2 and B = D, where D G is a symmetric positive definite 

matrix, the ;Bx{-eigenpairs are the D-eigenpairs defined in m, which satisfy 

Ax'^~^ = XDx, x^Dx = 1. (2-7) 
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• If -4. € m = m' and B = Z is the identity tensor, mode-1 ^S-eigenpairs are the 

eigenpairs defined in m. which satisfy 

( 2 . 8 ) 

where = [x^-\x^-\--- 

• If ^ € rI™’’’*], m = m' and B = X is the unit tensor, mode-1 ;Bij-eigenvalues are the 
H-eigenvalues defined in [37] . 

REMARK 2.3 Theoretical properties of mode-1 eigenvalues of tensors such as the Perron- 
Frobenius theory ( |21131 El 08] ) for nonnegative tensors can be parallelly developed for 
mode-fe eigenvalnes. However, as Horn and Johnson indicated in [18]: “One should not 
dismiss left eigenvectors as merely a parallel theoretical alternative to right eigenvectors. 
Each type of eigenvector can convey different information abont a matrix,” we believe 
that mode-1 through mode-m eigenpairs can convey different information about a general 
tensor of order m > 3. 


In the rest of this section, we will obtain an upper bound for the number of equivalence 
classes of mode-/c eigenpairs. As shown in Definition 12.11 Remark 12.11 and Remark 12.21 
the number of equivalence classes of mode-A; generalized eigenpairs for general tensors 
A € and B G Cl"*'’"] is equivalent to the nnmber of solutions to the following 

system of polynomials 




r(A,x) 


\aiXi -I- 02X2 -I-1- anXn + h) 


(2.9) 


where A and x := (xi,-- - ,x„)^ are the unknowns, ai,... ,an,b are random complex 
numbers. This motivates us to use Bernstein’s theorem and its extensions in the field 
of solving polynomial systems (see mm) to study the number of equivalence classes of 
eigenpairs. 

To initiate our discussion, we first introduce some commonly used notations and defi¬ 
nitions. Let P{x) := (pi(x),... ,pn{x))^ be a polynomial system with x := (xi,..., Xn)'^- 
For a := {ai ,..., an) G (Z^q)"^, write x" := x“^ • • • x"" and denote |q:| = ai -|- • • • J- an- 
Then P{x) can be denoted by 


P(x) 


/pi(x) := X] cp„x"\ 

oSSi 



E Cn,aX°‘ 
aeSn 


( 2 . 10 ) 


where Si,...,Sn are given finite subsets of (Z^q)"^ and Ci^a G C* := C\{0} are given 
coefficients of the corresponding monomials. Here for each i = 1,... ,n. Si is called the 
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support of Pi{x) and its convex hull Ri := conv(S'i) in M” is called the Newton polytope 
of Pi{x). {Si,..., Sn) is called the support of P{x). For nonnegative variables Ai,..., Xn, 
let Ai-Ri + • • • + XnRn be the Minkowski sum of XiRi, ■ ■ ■, XnRn-, be., 

XlRl + • • • + XnRn •= {Airi + • • • + XnVn \ Ti G Ri, i = 1, . . . , n}. 

The n-dimensional volume of Aii?i + • • • + XnRn, denoted by Vol„(Aiiii + ■ ■ ■ + XnRn), is a 
homogeneous polynomial function of degree n in Ai,..., A„ (See, for example. Proposition 
4.9 of [S] for a proof). The coefficient of the monomial A 1 A 2 ... A„ in Voln,(Aii?i + • • • + 
XnRn) is called the mixed volume of Ri,...,Rn, denoted by MXln{Ri, ■ ■ ■ ,Rn), or the 
mixed volume of the supports ^i,..., Sn, denoted by MV„(5i,..., Sn)- Sometimes it is 
also called the mixed volume of P{x) if no ambiguity exists. The following theorem relates 
the number of solutions of a polynomial system to its mixed volume. 

THEOREM 2.1 (Bernstein’s Theorem) [2] The number of isolated zeros in (C*)”, 
counting multiplicities, of a polynomial system P{x) = {pi{x),... ,Pn{x))'^ with supports 
Si,... ,Sn is bounded by the mixed volume MVn(S'i,..., Sn). Moreover, for generic choices 
of the coefficients in pi, the number of isolated zeros is exactly MV„,(S'i ,... ,Sn). 

An unexpected limitation of Theorem 12.11 is that it only counts the isolated zeros of a 
polynomial system in (C*)” rather than C”. To deal with this issue, Li and Wang gave 
the following theorem. 

THEOREM 2.2 [.lOj The number of isolated zeros in C”, counting multiplicities, of a 
polynomial system P{x) = {pi{x),... ,pn{x))'^ with supports Si,... ,Sn is bounded by the 
mixed volume MV„,(S'i U {0},... , S'n, U {0}). 

The following lemma was given as Exercise 7 on page 338 of [9]. 

LEMMA 2.1 Consider a polynomial system P{x) = {pi{x),... ,pn{x))'^ with supports 
Si = S 2 = ■■■ = Sn = S. Then 

MV„,(S ',... ,S) = n!Vol„(conv(S')). 

An upper bound for the number of equivalence classes of mode-A; eigenpairs which 
generalizes results in [31351 ED is given in the following theorem. 

THEOREM 2.3 Let A € and B G Assume that Bx"^' is not identically 

zero as a function of x. Let k be an integer such that 1 < k < m. Assume that A has 
finitely many equivalence classes of mode-k B-eigenpairs over C. 

(a) If m = m', then the number of equivalence classes of mode-k B-eigenpairs, counting 
multiplicities, is bounded by 

n{m — 1)”“^. 

If A and B are generic tensors, then A has exactly n{m — l)”"^ equivalence classes 
of mode-k B-eigenpairs. 
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(b) If m ^ m', then the number of equivalence classes of mode-k B-eigenpairs, counting 
multiplicities, is bounded by 

(m — I)"" — (m' — !)”■ 


m — m' 


If A and B are generic tensors, then A has exactly ((m — I)” — {m' — l)"’)/(m — m') 
equivalence classes of mode-k B-eigenpairs. 


Proof : Recall that the number of equivalence classes of mode-/c R-eigenpairs of A is equal 
to the number of solutions of (j2.9p . For the random hyperplane aixi + • • • + anXn + & = 0 
in (12.91) . without loss of generality, suppose that an ^ 0. Then Xn can be solved as 

Xn = CiXiI - V Cn-lXn-l + d, ( 2 - 11 ) 


where Cj = —aijan for z = — 1 and d = —bfan- Notice that the number of 

solutions of (12.9p in is the same as the number of solutions in C"" of the resulting 

system T*{X, xi,..., Xn-i) by substituting (j2.1ip into the first n equations of (12.9p . Denote 
the corresponding supports of T* by 5i,..., S'„. We claim that 


MV„(5iU{0},...,S„U{0}) < 


n(m — !)”■ 

(m — !)"■ — {m' — !)”■ 
m — m' 


m = m' 
m 7 ^ m' 


( 2 . 12 ) 


Let N denote the number of equivalence classes of mode-/c 0-eigenpairs of A over C. Then 
(12.121) implies that 


N < n{m — 1)"' ^ 


for m = m' and 

^ ^ (m — 1)*^ — (m' — I)"" 

“ m — m' 

for m ^ m'. When A and B are generic, equality holds in the two inequalities above by 
using Theorem 12.11 and Theorem 12.21 

To prove (I2.12p . let A € and B G be generic tensors. Similar to (|2.9I) the 

corresponding polynomial system to solve is 


/ (A(^)x^-^)i - A(Bx"^'-^)i \ 


f(X,x) 


{A^^^X^-^)n - XiBx^'-^)n 

\aiXi + 02X2 H-h OnXn + bj 


(2.13) 


Substituting (12.lip into the first n equations of (I2.13P yields a new system T*(A, xi,..., x„_i). 
The coefficients of the polynomials in the new system T* are the sums of products of cer¬ 
tain coefficients of the old system T. For example, when k = 1, the coefficient of in 
the first polynomial of the new system T* is 


m—1 


E 
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(2.14) 


















where is the set of all permutations of the set consisting of i numbers of n and m — 
1 — i numbers of 1. Since oi^o-’s and ci are all generic, the coefficient (I2.14p is nonzero. 
Similarly, all other coefficients in the new system are nonzero. Now let Si,... ,Sn be the 
corresponding supports of T*, we can assume that all monomials 


rr°‘lrr°‘2 


n 


Oii G ^>Q, Oil CH2 “b * * * “1“ Q^n — m — 1} 


and 




“1 

2 


ai G Z>o, ai + 02 H- h an = m' - 1} 


will appear in each of the hrst n equations in (I2.13p . Therefore, after substituting (12.lip 
into the first n equations of (I2.13h . all monomials 




0(2 


. . . X 


0!n-l 

n—1 


ai G Z>o, ai + q ;2 H-b an-i < m - 1} 


and 

{Ax“^X 2 ^ ... ai G Z>o, ai + 0:2 H-b On-i < m' - 1} 

will be contained in each equation of T*. This implies that Si,... ,Sn are all equal to 


S := {(0, a) I a G (Z>q^)^, |a| < m — 1} U {(1, a) | a G (Z>q^)^, |a| <m! — 1}. 

Notice that the convex hull of the set {a G (Z”g^)^ | |a| < m — 1} is the (n — l)-simplex 
in with vertices (0,0, • • • , 0), (m — 1,0, • • • , 0), • • •, (0, • • • , 0, m — 1), and the con¬ 
vex hull of the set {a G (Z”g^)^ | |q;| < m' — 1} is the (n — l)-simplex in with 

vertices (0,0, ,0),(m' — 1,0,--- ,0),--- , (0, • • • ,0,m' — 1). Thus, their volumes are 

(m — l)"'“^/(n — 1)! and {m' — l)"'“^/(n — 1)! respectively (see, for example. Exercises 2 
and 3 on page 307 of [9]). Let Q be the convex hull of S. Then Q is the linear interpolation 
between the two aforementioned simplices with 


Q = {(t, conv(a)) 11 G [0,1], a G |q;| < m — 1 -b t{m' — m)}. 

Since for fixed t, conv({a G (Z^g^)"^ | |a| < m — 1 -b t{m' — m)}) is a simplex with volume 
(m — 1 -b t{m' — m))"'“^/(n — 1 )!, we have 


Voln 



(m 


1 -b {m' — m)t)"' ^ 
(n — 1 )! 


dt 


' {m — !)"■ — (m' — 1 )*^ 

(m — m')n\ 
n{m — 1 )”“^ 


m 7 b rn', 
m = m'. 


Therefore, by Lemma l2.ll 


MV„(5i,...,5„) = n!Vol„(g) 


(m — 1 )” — {m! — !)"■ 

m — m' 
n{m — l)”"”^. 


m 7 b rn', 
m = m'. 


When tensors A and B in (|2.9p have some zero entries, each support Si of T* must 
be a subset of the support Si of T*, i.e.. Si C Si. In addition, since each polynomial of 
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T* contains a constant term (which arises when substituting (I2.1ip into the term in 

the first n polynomials of T), 0 G 5*. Thus we have Si U {0} C Si. Since mixed volume is 
monotonic [1], 

MV„(5i U {0},..., U {0}) < MV„(5i,..., 5„). 


This implies that (j2.12D holds. 


□ 


REMARK 2.4 A few remarks about Theorem 12.31 

(a) Theorem 12.31 provides an upper bound for the number of equivalence classes of B- 
eigenpairs of tensor A if A has finitely many B-eigenpairs. The bound is tight when 
A and B are generic. In the literature, the numbers of eigenvalues as defined in (j2.8p 
and E-eigenpairs have been investigated. 

When m = m' and B is the mth order, n-dimensional identity tensor, Qi [371 Theo¬ 
rem 1] and Chang, Qi, and Zhang [U Remarks 1] proved that a tensor A € has 

n(m—1)”“^ eigenvalues (12.8p . In this case, Part (a) of Theorem [2A] gives re(m —1)"'“^ 
as the upper bound for the number of equivalence classes of such eigenpairs. 

When m' = 2 and B is the n x n identity matrix, Cartwright and Sturmfels [U 
Theorem 1.2] proves that tensor A has ((m — 1)"' — l)/(m — 2) equivalent classes of 
E-eigenpairs if it has finitely many E-eigenpairs. Part (b) of our Theorem 12.31 gives 
((m — 1)"^ — l)/(m — 2) as the upper bound of the number of equivalence classes of 
E-eigenpairs if A has finitely many E-eigenpairs. 

(b) The upper bound given in Theorem 12.31 can be highly useful in designing effective 
homotopy methods for computing mode-fc generalized eigenpairs. In fact, the ho- 
motopy method described in Algorithm 13. II for the case m = m' relies on the bound 
n(m — 1)”“^. 

3 Computing complex tensor eigenpairs via homotopy meth¬ 
ods 

Consider A G and B G As discussed in Section [2l the problem of computing 

mode-/c R-eigenpairs of A in (|2.3jl or (12.411 is equivalent to the problem of solving (12.91) . 
and if m A normalize (A,x) to satisfy that Rx™' = 1. Since ()2.9I) is a polynomial 
system, we consider to use a homotopy continuation method to numerically solve it. 

The basic idea of using homotopy continuation method to solve a general polynomial 
system P{x) = (pi(x),... ,pn{x))'^ = 0 as defined in (I2.10p is to first deform P{x) = 0 
to another polynomial system Q(x) = 0 that is easy to solve. Specihcally, we construct a 
homotopy H -.PPx [0,1] —>■ C"" such that H{x, 0) = Q{x) and H{x, 1) = P{x). Then under 
certain conditions, the homotopy H{x,t) = 0 has smooth solution paths parameterized by 
t for t G [0,1) and all the isolated solutions of P{x) = 0 can be reached by tracing these 
paths. 
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A useful homotopy is the linear homotopy (see [Ml EH sa SS]): 

H{x, t) = (1 - t)jQ{x) + tP{x) = 0, t e [0,1], 


(3.1) 


where 7 is a generic nonzero complex number. It is very critical to choose a suitable Q{x) 
such that the system Q{x) = 0 is easy to solve and all isolated solutions of P{x) = 0 can 
be found by tracing solution curves of H{x,t) = 0. 

One choice of Q{x) that always makes the linear homotopy ()3.ip work is the so-called 
total degree homotopy, in which the starting system Q{x) = 0 has deg = di x d 2 x ■ ■ ■ x dn 
solutions ([Ml EH 05]), where di,... ,dn are the degrees of polynomials pi{x),... ,Pnix) 
respectively, deg is called the total degree or Bezout’s number. By tracking deg number 
of solution paths of ()3.ip we can find all the isolated solutions of P{x) = 0. However, 
most polynomial systems in applications usually have far fewer than deg solutions. In this 
case, many of the deg paths will diverge to infinity as t —>■ I resulting in huge wasteful 
computations. 

The polyhedral homotopy [22] based on Bernstein’s Theorem [2] makes significant 
progress in this sense. In this method, the number of paths that need to be traced is the 
mixed volume of a polynomial system, which generally provides a much tighter bound 
than Bezout’s number for the number of isolated zeros of a polynomial system. Hence 
the new method reduces a significant amount of extraneous paths than the total degree 
homotopy in most occasions and thereby is much more efficient. However, the polyhedral 
homotopy includes two major stages: mixed volume computation and tracking paths. The 
computation of mixed volumes is a sophisticated procedure [3|. Moreover, mixed volume 
computation can be very expensive for large polynomial systems. Thus if the mixed 
volume is far less than the Bezout’s number and an appropriate linear homotopy can be 
constructed so that only mixed volume number of paths need to be traced, the system is 
better to be solved by using a linear homotopy instead of the polyhedral homotopy. 

To compute tensor eigenpairs, one can certainly use the polyhedral homotopy imple¬ 
mented in HOM4PS [M], PHCpack [S], PHoM [H], PSOLVE [39] (which is a MATLAB 
implementation of HOM4PS), or the total degree homotopy implemented in Bertini [3]. 
However, using these methods to solve (|2.9I) or (12.41) directly does not take advantage of 
the special structures of a tensor eigenproblem. We will introduce two homotopy-type 
algorithms here that utilize such structures. 

3.1 A linear homotopy method when m = m! 

Theorem 12.31 gives us that the mixed volume of (j2.9p when m = m' is n{m — I)”~^, which 
is far less than the Bezout’s number, m”". We consider constructing a linear homotopy in 
which only the mixed volume number of paths are traced. 

For a polynomial system P{x) = {pi{x),... ,pn{x))'^ as defined in (|2.I0p . where x = 
(xi,..., Xn). Partition the variables xi,... ,Xn into k groups yi = (x^^\ ..., x[^^),y 2 = 

(x^^^, • • •, x\‘^'^),... ,yk = , ■ ■ ■, x\^'^) with li + ■ ■ ■ + Ik = n. Let dij be the degree of pi 

with respect to yj for i = I,...,n and j = 1,... ,k. Then the multihomogeneous Bezout’s 
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number of P{x) with respect to (yi,..., is the coefficient of ™ product 

n 

Y[{ii lai + • • • + dikCtk)- 
i=l 

The following theorem will play an important role in constructing a proper linear 
homotopy. 


THEOREM 3.1 [l3] Let Q{x) be a system of polynomials chosen to have the same mul- 
tihomogeneous form as P{x) with respect to certain partition of the variables (xi,..., x„). 
Assume Q(x) = 0 has exactly the multihomogeneous Bezout’s number M of nonsingular 
solutions with respect to this partition. Then for almost all 7 G C*, the homotopy 

H{x,t) = (1 - t)^Q{x) +tP{x) = 0, 

has M nonsingular solution paths on t ^ [ 0 , 1 ) whose endpoints as t —)■ 1 include all the 
isolated solutions of P{x) = 0. 


For (j2.9p . when m = m' the following polynomial system 


/ - \{Bx^-% \ 


G{\x) 


\aiXi + 02X2 H-+ OnXn + b) 


needs to be solved, where A and x := (xi, • • • ,x„)^ are the unknowns, oi, 
random complex numbers. Consider the starting system 


Q(A,x) 


/(A-/ii)(xr'-/?i)\ 

(A-/i2)(a;r'-^2) 

(A-/r„)(xr'-/3n) 

\ CiXi + . . . CnXn + d ) 


(3.2) 


On,b are 


(3.3) 


where Hi, fdi, Ci for i = 1,... ,n and d are random nonzero complex numbers. 


THEOREM 3.2 Let G{\,x) and Q{\,x) be defined as 113.^) and il3.3\) respectively. Then 
all the isolated zeros (A,x) in of G{\,x) can be found by using the homotopy 

H{\,x,t) = (1 — t)^Q{\,x) + tG{X,x) = t), f G [0,1] (3-4) 

for almost all 7 G C*. 

Proof : It is sufficient to verify that Q{\,x) satisfies all the assumptions of Theorem 13.11 
Partition the variables (A, xi,..., Xn) into two groups: (A) and (xi,..., Xn), we can easily 
see that each of the first polynomial equations in (j3.2p and (13.3h has degree 1 in (A) and 
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degree m — 1 in (xi,..., Xn), and the last equation in both (13.211 and (I3.3p has degree 0 in 
(A) and degree 1 in (xi,... ,Xn)- Hence (|3.2I1 and (13.3p have the same multihomogeneous 
Bezout’s number as the coefficient of the product 

[1 • ai + (m - l)a2]”(0 • ai + 1 • 02 )- 

It can be easily computed that this coefficient is equal to 

(m- =nim-l)^-\ 

Hence ()3.2I1 and (|3.3p have the same multihomogeneous Bezout’s number n{m — 1)”“^ 
with respect to the partition (A) and (xi,..., x„). 

We now show that Q{X,x) in (13.311 has exactly n{m — 1)”“^ zeros. Notice that if A 
is equal to none of fii,... ,^n, then we end up with a system of n + 1 equations and n 
unknowns, which has no solutions. Thus A must be equal to one of //!,...,//„• Without 
loss of generality, assume that A = /ii. Then xi,..., x^ can be determined by 

x^n-i _ p, ^ i = 2 ,...,n 

CiXi + . . . CnXn + d = 0 

Obviously, each x* for i = 2,... ,n can be chosen as one of the (m — l)-th root of /?* and 
xi will be solved by substituting the chosen X 2 ,..., x^ into the last hyperplane equation. 
So there are (m — 1)"’“^ solutions corresponding to A = /ii. This argument holds for A 
being any //*. Therefore, there are totally n{m — solutions. 

It remains to prove that each solution of Q{\,x) = 0 in (|3.3I1 is nonsingular. As 
discussed above, any solution (A*,x*) of (j3.3h satisfies 

A* = //i, 

{x*)^~^-Pj = 0, j = I,--- ,i-l,i + l,--- ,n, (3.5) 

c\Xi + • • • + CnX.^ d = 0. 

Let DQ{X,x) be the Jacobian of Q{X,x) with respect to (A, x). It is sufficient to show 
that DQ{X*, X*) is nonsingular. Denote 

Bj{X,x) := (A - - l)x™“^ 

\ 

Bi-i 

Bi 

Bi+i 


Q —1 Ci Cj+l . . . Cn j 


Aj(X,x) 

for j = 1,... ,n. Then 


DQiX,x) = 




/ Ai Hi 

Aj-i 

Ai 

Aj+i 

Aji 

Vo Cl 
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Note that 74j(A*,x*) = — f3j =0, j i and Bi{X*,x*) = (A* — Hi){m — 

l)(x*)™'“^ = (/ij — — l){x*)^~^ = 0 by p.Sp . For simplicity, write Aj := Aj{X*,x*) 

and Bj := Bj{X*,x*). Then 

/O Bl \ 


DQ{X*,x*) 


0 

0 


B 


* 

2—1 


0 


B 


* 

i+1 


Then 


0 B*^ 

\ 0 Cl ... C^—i Ci CiJ^i ... Cfi J 

det{DQ{X*,x*)) = i-iy+^A*i-l)^+^CillB* / 0 


by (1331) . 


□ 


Theorem 13.21 suggests us that (j3.4p can be used to solve ^2.90 in the case of m = m'. 
For simplicity, write u := (A,x). In order to improve numerical stability, we apply the 
transformation s = Inf to (|3.4I) (a strategy first suggested in [T3]) and obtain the new 
homotopy as 

H{u, s) = (1 — e^)'yQ{u) + e^G{u) = 0, s € [—oo, 0] (3-6) 

where Q{u) = Q{X,x) and G{u) = G{X,x) are defined in (|3.3p and (|3.2I) respectively. 

We now introduce our linear homotopy method for computing mode-A: generalized 
eigenpairs when m = m!. 

ALGORITHM 3.1 (Compute complex mode-A: ;B-eigenpairs of A, where A,B & Cl™’"'!.) 


Step 1 . Compute all solutions of Q{u) as defined in I13.3\} . 

Step 2. Path following: Follow the paths from s = —oo to s = 0. In reality, we 
certainly cannot start from s = —oo. In this case, one can choose a very negative sq <ind 
obtain a starting point by using Newton’s iterations: 

-[Hu{w^’^\so)]-^II so), A: = 0,1,... 


until so)|| is very small for some N. Here is a solution of Q{u) = 0. Let 

no := n(so) and take uq = . Then path following can be triggered. 

Path following is done using the prediction-correction method. Let {uk, Sk) '■= {u{sk), Sk), 
to find the next point on the path H{u, s) = 0, we employ the following strategy: 


du — 

Prediction Step: Compute the tangent vector — to H{u,s) 


linear system 


IIu{,Uk,Sk) , — Hgl^Uk, Sk) 

ds 


0 at Sk by solving the 
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. du _ , 

for Ihen compute the approximation u to u^+i by 
ds 


u = Uk + As 


du 
ds ’ 


•Sfc+i — Sk + As, 


where As is a stepsize. 


• Correction Step: Use Newton’s iterations. Initialize vq = u. For i = 0,1,2,..., 
compute 

Vi+i =Vi- [Hu{vi,Sk+i)]~^H{vi,Sk+i) 
until sj)|| is very small. Then let Uk+i = vj. 

Step 3. End game. During Step 3 when sn is very close to 0, the corresponding um 
should be very close to a zero u* of G{u) = G{X,x). So Newton’s iterations 

^{k+i) ^ ^(k) _ A: = 0,1,... 

will be employed to refine our final approximation u to u*. If DG{u*) is nonsingular, then 
u will be a very good approximation of u* with multiplieity 1. If DG{u*) is singular, u 
is either an isolated singular zero of G{u) with some multiplieity I > 1 or in a positive 
dimensional solution component of G{u) = 0. We use a strategy provided in Chapter VIII 
of J2^ (see also f43\l ) to verify whether u is an isolated zero with multiplicity I > 1 or in 
a positive dimensional solution eomponent of G{u) =0. 

Step 4. For eaeh solution u = (A, x) obtained in Step 3, normalize x in the following 
fashion to get a new eigenveetor 

y = ^ (3-7) 


^20 


can be obtained, where i^ := arg|xj|. Hence (A,y) is an eigenpair. 

REMARK 3.1 A few remarks about Algorithm l3.lt 

(a) As defined in (|2.3|) and Remark 1 2 . It if (A,x) is an eigenpair, {\,tx) for f 7 ^ 0 is also 
an eigenpair. Therefore, Step 4 is well defined in this sense. 

(b) Notice that if x is a real eigenvector associated with a real eigenvalue A, tx for any 
t E C\{0} will be a complex eigenvector associated with the same eigenvalue A. If 
in any case a complex eigenvector like tx is obtained in Step 3 of Algorithm 13.11 
applying (13.71) to tx will give us a new real eigenvector. In this sense. Step 4 is very 
helpful for us to detect real eigenpairs. 

(c) According to Theorem 12.3t if A has finitely many equivalence classes of R-eigenpairs, 
then the number of equivalence classes of R-eigenpairs, counting multiplicities, is 
bounded by n(m—1)”“^. Moreover, this bound is attained when A and B are generic. 
This result implies that the optimal number of paths to follow in a homotopy method 
for solving the system (|3.2I) is n{m — l)””^. Our starting system (13.31) has exactly 
n{m — 1)"'“^ nonsingular solutions. In this sense. Algorithm 13.11 follows the optimal 
number of paths for solving the system (13.21) . 
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3.2 A polyhedral homotopy method when m ^ m! 

To compute mode-fe generalized tensor eigenpairs when m ^ m', we use the equivalence 
class structure of the eigenproblem as described in Remark [2.II We first solve (12.91) to hnd 
a representative (A, x) from each equivalence class and then find all m' eigenpairs from 
each equivalence class by simply using \' = \^x' = tx, where t is a root of t"* = 1 . 

We use a polyhedral homotopy method to solve the system ()2.9p . In our implementation, 
the polyhedral homotopy method is PSOLVE ([l9]), with some modihcations, as described 
in Subsection 3.3. The modifications are based on Strategies 2 and 3 introduced in the 
next subsection. 

One may think of solving ( 12 .4p directly to get m' eigenpairs from each equivalence 
class. However, this alternative method would have to follow m! times as many paths as 
the approach we described in the previous paragraph and therefore it would need much 
more computation. 

Now we present our algorithm for computing mode-A: generalized eigenpairs when m ^ 
m!. 

ALGORITHM 3.2 (Compute complex mode-A; H-eigenpairs of A, where A € , B € 

with m A m'.) 

Step 1. Using modified PSOLVE to get all solutions (A,x) of 112.9\) . 

Step 2. For each (A, x) obtained in Step 1, if Hx™ 7 ^ 0, normalize it to get an 
eigenpair (A*,x*) by 


to satisfy 12.4\ ). 

Step 3. Compute m' equivalent eigenpairs (A',x') of {X*,x*) by X' = A* and 

x' = tx* with t being a root of = 1. 

3.3 Implementation tips 

By using random complex numbers in the formulation of homotopy, theoretically, with 
probability one the solution paths do not cross or go to infinity in the middle. In practice, 
however, two paths may become very close to each other and the magnitude of some com¬ 
ponents of a solution curve may become very large during the procedure of path tracking. 
This causes various numerical difficulties such as missing solutions or losing efficiency or 
stability. In our implementation of Algorithms 13. II and 13.21 we use the following strategies 
to address these issues. We focus our discussion on Algorithms 13.11 Step 1 of Algorithm 
M\ is done similarly. 

When tracing two paths that are sufficiently close, it is possible for the path tracing 
algorithm to jump from one path to the other path and thus result in the missing of zeros. 
To minimize the chance for curve jumping and keep the efficiency, our First Strategy 
is: The stepsize As in Step 2 of Algorithm \S.1\ is chosen adaptively. Initially, set sq = 
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—20(n + 1), where n + 1 is the number of unknown variables A, xi,..., in (|2.9p or (13. 2p . 
and As = — so/3. Similar to [25], if more than 3 steps of Newton iterations were required 
to converge within the desired accuracy, As is halved and the shorter step is attempted. 
On the other hand, if several (the default being 2) consecutive steps were not cut. As is 
doubled, up to a prescribed maximum value (the default being — Sfc/3). 

Although this adaptive approach can often significantly reduce the possibility of curve 
jumping, it can still occur in some cases. Our Second Strategy is: To check if there is 
curve jumping, we store all the found solutions in a binary search tree. Each time when a 
new solution is found, we can quickly find (with time complexity OilogN), where N is the 
number of solutions) whether there is any existing solution that is numerically identical 
to the new solution, that is, the difference between them is less than a threshold (the 
default being If two numerically identical solutions are detected and the condition 

numbers of their Jacobian matrices are greater than a threshold (the default being 10^® 
we consider that curve jumping likely has occurred. We then retrace the two associated 
curves with more restrictively chosen parameters in the projective space, as described in 
the paragraph after the next one. 

When the magnitude of some components of certain solution curves become very large 
in the middle, tracing these paths may fail due to numerical instability. This issue can 
be largely resolved by following paths in the projective space (see, for example, [33]). 
However, empirically it is more time consuming to follow all paths in the projective space 
than in the complex space. In our implementation of Algorithms 13.11 and 13.21 our Third 
Strategy is: To retrace solution curves in the projective space only for those paths that are 
detected to have very large solution components. 

To trace a path in the projective space, we first homogenize each polynomial equation 
of (13. 6 |) in the variables A, xi,..., x„ to get the homotopy 





/ xo(A'‘A’^-fi-X(Bx"‘-^)i \ 

H{X,x,s) = (1 - 

(A - fJ,nXo){x'J~^ - JnX^~^) 

\ ClXl + . . . CnX„ + dxo / 

1 s 

+ e 

xo(A^A'^-^)u - X(13x’"-^)„ 
VaiXi + 02X2 H-h a„x„ + bxo/ 


(3.8) 

where x = (xo,xi,... ,x„)^, and then follow the solution curve of (13.Sp in the projective 
space. Notice that in (13.8p if (A, x) is a solution, so is (aA, ax) for a € C\{0}. Thus along 
the path we can always scale (A,x) to keep each component’s magnitude in a suitable 
finite range. 

To the best of our knowledge. Strategies 2 and 3 have not been used in other imple¬ 
mentations of homotopy methods, although some packages may trace all curves in the 
projective space. 

4 Computing real tensor eigenpairs via homotopy methods 

In some applications, tensor A is real and only real eigenpairs (or real eigenvalues) of 
A are of interest ( |12l 137] ). In this situation, only real zeros of the polynomial system 
(12.4p or (|3.2h are needed. It is worth noting that there is currently no effective method 
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to find all real zeros for a polynomial system directly. One may think to use a homotopy 
continuation method to trace only real zeros from the start system to the target system. 
However, this approach generally does not guarantee a real zero at the end, because the 
homotopy methods inherently have to trace paths in the complex space in order to avoid 
the discriminant locus. 

For a real tensor A, a real eigenvalue may have complex eigenvectors. Sometimes 
identifying which eigenvalues are real is the only concern. In this case, we can first compute 
complex zeros (A,x) of ()3.2p by Algorithm 13. 1 1 or (j2.4l) by Algorithm 13.21 then identify the 
real eigenvalues by checking the size of the imaginary parts of A’s. Specifically, let (A*, x*) 
be a computed eigenpair. If 

|Im(A*)| <<5o, 

where Sq is a threshold for the imaginary part (the default value 6o = 10 “®), then we take 
Re(A*) as a real eigenvalue. 

Note that if is a nonzero integer multiple of 4 (for example, m = 5, = 4 or 

m = 10 , m' = 8,) and A has an eigenpair (A*, x*) with a purely imaginary eigenvalue A* = 
bi, where 6 G M, then one can easily show that ( 6 , (—and (—5, 
are eigenpairs with real eigenvalues. Therefore, when is a nonzero integer multiple 

of 4, if (A*,x*) is an eigenpair found by Algorithm 13.21 such that 

|Re(A*)| <<5o, 

then we take Im(A*) and — Im(A*) as real eigenvalues, with corresponding eigenvectors 
and 

When looking for real tensor eigenpairs (i.e., both eigenvalues and eigenvectors being 
real), the situation becomes more complicated. We use a two-step procedure. The first 
step is to get complex zeros (A, x) of (13.2p by Algorithm 13.11 or (12.4p by Algorithm 13.21 
The second step is to extract all real eigenpairs (A, x) from the complex zeros obtained in 
the first step. 

To facilitate the discussion, the following notation is introduced. For a vector a = 
(ai,... jOn)^ G C, let 

Im(a) = (Im(ai),..., Im(a„))'^, Re(a) = (Re(ai),..., Re(a„))^. 

Suppose that (A*,x*) is an eigenpair found in the first step. Consider two cases: 
(i) (A*,x*) is an isolated eigenpair; (ii) (A*,x*) is an eigenpair contained in a positive 
dimensional solution component of system (|3.2I) or (12.4p . 

When (A*,x*) is an isolated eigenpair, take (Re(A*), Re(x*)) as a real eigenpair if 

II Im(A*,x *)||2 < 5o- 

If (A*, X*) is an eigenpair in a positive dimensional solution component of system (j3.2p 
or (|2.4I) . in general real eigenvectors are not guaranteed to be found by Algorithm 13.11 
or Algorithm 13.21 even if the corresponding eigenvalue A* is real. In this case, we will 
construct the following Newton homotopy [T] 

i/(A,x,t) := P(A,x) - (1-t)P(A*,Re(x*)), t G [0,1] (4.1) 
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to follow curves in (A, x) € in order to get a real eigenpair. Notice that when following 
curves in the complex space it is proved in [2^ that the solution curves of ()4.1I) can be 
parameterized by t, but the solution curves of (|4.1|) is not necessarily to be a function of t 
when restricted in the real space. So a different method to follow curves is needed. In this 
case parameterizing the solution curves by the arc length s is suggested in m- Instead 
of following paths using (|4.ip . we will use the homotopy 

H{{X{s),x{s),t{s)) =0. (4.2) 

For a description of the Newton homotopy method, we refer to m- 

An interesting phenomenon we have observed in our experiments is that in some 
cases, the real eigenpairs can be obtained more straightforwardly from the complex eigen- 
pairs found from Algorithm 13.11 or Algorithm 13.21 If (X*,x*) is in a positive dimen¬ 
sional solution component of (|2.4I) and A* G R, then (A*,Re(x*)/(. 8 Re(x*))^/™’^) and 
{X*,lm{x*)/{Blm{x*))^/^') can be mode-/c Br eigenpairs of A. This gives us a heuristic 
approach to find real eigenpairs for eigenpairs belong to positive dimensional components. 
We remark that this approach works well for all the examples (e.g.. Example 4.8, 4.11, 
4.13, 4.14) in [12] when a real Z-eigenvalue has infinitely many real Z-eigenvectors. The 
following Proposition gives a justification for the approach in special cases. 

PROPOSITION 4.1 Let A G r[™’" 1 and B G R[™'^’"'1. Let k be an integer such that 
\ < k < m. Let X G'R he a real mode-k B eigenvalue of A. If U := {x G C"' | = 

XBx^ contains a complex linear subspace V of such that y G V implies y G V , 
then for any x = f, + ip G V such that r/ G R” and / 0, r/ / 0, 

• When m = m', ^ and rj are both real mode-k B eigenvectors of A associated with A. 

• When m A 7^ 0 and Brf^ A 0; ^he normalized vectors 

?; - in '= - 

are real mode-k B eigenvectors of A associated with A. 


Proof : Let x G V. Then x G V. Since R is a linear subspace, ^ = (x -|- x)/2 and 
7] = (x — x)/{2i) are also in V. Thus, when m = m', ^ and r] are both real mode-/c B 
eigenvectors of A associated with A. If m 7 ^ m', we have 


Bv^' = Bi^i2-i^WhVi2 ■ ■ ■ 


E « 








B^m' 


BC 


B^m' 


= 1 . 
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This implies that u is a real mode-fc B eigenvector of A associated with A. Similarly, 
Bw^ = 1 can also be verified. Therefore, w is also a real mode-A: B eigenvector of A 
associated with A. □ 


REMARK 4.1 A natural question is when U dehned in Proposition 14. II contains a linear 
subspace V. Consider the case when A is a symmetric tensor with a low rank decompo¬ 
sition. For simplicity, consider m = 3. For vectors a,b,c € C”", define the outer product 
tensor aoboc = {uibjCk)- Suppose that a symmetric tensor A G Ml™’”! can be decomposed 
as 

= yi o yi o yi 7/2 o 2/2 o 2/2 H-h 2/r o 2/r o 2/r, 

where r < n, t/k & K"', k = 1,2, ■ ■ ■ ,r. 

Let W = span( 2 /i, 2 / 2 , • • • > Ur) and let V be the orthogonal complement of W. Then V 
is a linear subspace of C”'. Clearly, x G F implies x G F. Moreover, for any x G F\{0}, 

Ax'"-^ = 0. 

Thus X is an eigenvector of A corresponding to the eigenvalue 0. In this case, Re(x)/|| Re(x) || 
and Im(x)/|| Im(x)|| are real Z-eigenvectors of A corresponding to the real Z-eigenvalue 0. 
Hence, the set F is a desired linear subspace of C"" contained in U. 

Finally we present an algorithm for computing real eigenpairs that combines the heuris¬ 
tic approach and the Newton homotopy method. 

ALGORITHM 4.1 (Compnte real mode-/c R-eigenpairs of A, where A G G 

kK.A) 


Step 1 . Compute all complex eigenpairs using Algorithm \3.1\ or Algorithm \3.A Let 
K denote the set of found eigenpairs (A, x) sueh that |Im(A)| < Sq. 

Step 2. For eaeh eigenpair (A*,x*) G K: if (A*,x*) is in a positive dimensional 
solution eomponent of IS. gj) or \2.4^ , go to Step 3. Otherwise, (A*,x*) is an isolated 
eigenpair. If || Im(x *)||2 < Sq, then take (Re(A*),Re(x*)) as a real eigenpair and stop. 

Step 3. Set A = Re(A*). If m = m', set v := Re(x*) (ifKe{x*) A Oj and w := Im(x*) 
(iflm{x*) 7 ^ 0); otherwise, set 


Re(x*) 

(RRe(x*)”^')Vm' 


(if RRe(x*)™V0), 


and 


Im(x*) 

(RIm(x*)”^')V"*' 


(if Blm{x*r' AO). 


If (A, v) or (A, w) is a mode-k B-eigenpair of A, then we have obtained a real eigenpair 
and stop. Otherwise, goto Step 4- 

Step 4. Starting from (A*,x*), use the Newton homotopy method to follow curves of 
i4.^ to find a real eigenpair. 
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5 Numerical results 


Based on the algorithms introduced in Section 3 and Section 4, a MATLAB package 
TenEig has been developed. The current version is TenEig 1.1. The numerical results 
reported in this paper were obtained using this version. The package can be downloaded 
from 


http://www.math.msu.edu/~chenlipi/TenEig.html 


Consider the tensors A G and 13 G In the TenEig package, function teig 

can be used to compute the general mode-A: B eigenvalues and eigenvectors of a tensor 
A for m = m!. The input of this function is: tensor A or the polynomial form (if A is 
symmetric) Ax^, tensor B (the default is the identity tensor), mode k (the default value 
is 1), and the output is: mode-Zc B eigenvalues and eigenvectors of A. By default, teig 
hnds eigenvalues and eigenvectors for the eigenproblem (j2.8p . 

The function teneig computes the general mode-Zc B eigenvalues and eigenvectors of 
a tensor A for m A The input of this function is: tensor A or the polynomial form 
(if A is symmetric) Ax"^, tensor B, mode k (the default value is 1), and the output is: 
mode-A: B eigenvalues and eigenvectors of A. If B is chosen as the identify matrix, the 
teneig computes the E-eigenvalues and E-eigenvectors of A as defined in Qi m- 

Since E-eigenpairs of a tensor are frequently needed, our package includes a separate 
function eeig, which only computes E-eigenpairs of a tensor. 

The package also includes two functions heig and zeig to compute real eigenpairs of 
a tensor: The first one computes H-eigenpairs and the second one computes Z-eigenpairs. 

In the next two subsections, numerical results are reported to illustrate the effectiveness 
and efficiency of our methods for computing tensor eigenpairs. All the numerical experi¬ 
ments were done on a Thinkpad T400 Laptop with an Intel(R) dual core CPU at 2.80GHz 
and 2GB of RAM, running on a Windows 7 operating system. The package TenEig was 
run using MATLAB 2013a. In our examples, we used teig or teneig to compute gener¬ 
alized eigenpairs, teig to compute eigenpairs (j2.8j) . eeig to compute E-eigenvalues, heig 
to compute (real) H-eigenpairs, and zeig to compute (real) Z-eigenpairs, respectively. 

5.1 Examples for computing complex eigenpairs 

In this subsection, some numerical examples illustrating the performance of TenEig for 
computing complex tensor R-eigenpairs are provided. 

A numerical solver NSolve, based on the Grobner basis, for solving systems of algebraic 
equations is provided by Mathematica. We will compare the performance of TenEig and 
NSolve in computing complex tensor H-eigenpairs. Denote 

T{m,n) := n(m —1)”“^, 

E{m,n) := ((m — I)"" — l)/(m — 2), 

G{m,m',n) := ((m — I)"" — (m'— l)"')/(m — m'). 
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Recall that Theorem 12.31 explains that for tensors A € and B € the number 

of equivalence classes of isolated 0-eigenpairs of A is bounded by T(m, n) for m = m' and 
G{m,m',n) for m A particular, Remark 12.41 states that the number of equivalence 

classes of isolated eigenpairs of the eigenproblem (j2.8p is bounded by T{m,n) and the 
number of equivalence classes of isolated E-eigenpairs of A is bounded by E{m,n). 


EXAMPLE 5.1 In this example, we compare the performance of our TenEig with 
NSolve and PSDLVE. teig, eeig, NSolve and PSDLVE are used to compute eigenpairs 
(12.8p or E-eigenpairs of a generic tensor A G We remark the following: 

(a) teig is based on Algorithm 13.11 The polynomial system solved is (|3.2p . 

(b) eeig is based on Algorithm 13.21 The polynomial system solved is (j2.9p . 

(c) When computing eigenpairs (12.81) . NSolve and PSDLVE solve the polynomial system 
defined by (13.2|) . 

(d) When computing E-eigenpairs, NSolve and PSDLVE solve the polynomial system 
defined by (12.4p . 

The tensors A were generated using randn{n, • • • , n)+i*randn(n, ■ ■ ■ , n) in MATLAB. 
The results are given in Tabled) In this table, N denotes the number of equivalence classes 
of eigenpairs (12.81) or E-eigenpairs found by teig, eeig, PSDLVE or NSolve, the reported 
CPU times are in seconds, denotes that no results were returned after 12 hours. 


(m, n) 

T{m, n) 

Alg 

N 

time (s) 

E{m, n) 

Alg 

N 

time (s) 



teig 

405 

15.8 


eeig 

121 

5.4 

(4,5) 

405 

PSDLVE 

404 

14.0 

121 

PSDLVE 

121 

9.5 



NSolve 

405 

3136.4 


NSolve 

121 

486.6 



teig 

1280 

73.8 


eeig 

341 

22.3 

(5,5) 

1280 

PSDLVE 

1280 

65.5 

341 

PSDLVE 

341 

38.6 



NSolve 

- 

- 


NSolve 

341 

9264.8 



teig 

6144 

606.5 


eeig 

1365 

166.5 

(5,6) 

6144 

PSDLVE 

6144 

694.2 

1365 

PSDLVE 

1365 

283.6 



NSolve 

- 

- 


NSolve 

- 

- 



teig 

18750 

3721.3 


eeig 

3906 

990.2 

(6,6) 

18750 

PSDLVE 

18748 

4636.0 

3906 

PSDLVE 

3905 

1721.0 



NSolve 

- 

- 


NSolve 

- 

- 


Table 1: Comparison of teig and eeig with PSDLVE and NSolve 

Erom table dl we see that teig and eeig successfully find all equivalence classes of 
eigenpairs (12.81) or E-eigenpairs using reasonable amount of time in all the cases. NSolve 
cannot get any results in 12 hours in some cases (we terminated it after 12 hours). Al¬ 
though PSDLVE successfully finds all equivalence classes in many cases, it does miss a few 
equivalence classes in some cases. We believe that the robustness of teig and eeig is due 
to their use of the retracing strategies and working in the projective space, as described 
in Subsection 3.3. Regarding the CPU time usage, PSDLVE is comparable to teig when 
the number of eigenpairs ()2.8[) T(m, n) is moderate. When T{m,n) gets larger, teig uses 
less time. This is because the mixed volume computation in PSDLVE takes significantly 
more time when T(m, n) becomes large. We also observe that eeig uses less time than 
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PSOLVE. Note that by using the equivalent class structure discussed in Remark l2.ll eeig 
traces E{m,n) paths. On the other hand, PSOLVE traces 2E{m,n) paths because here it 
solves the system (12.41) from Definition 12.11 directly. 

EXAMPLE 5.2 In this example we show the effectiveness and efficiency of teig for 
finding all equivalence classes of isolated eigenpairs as defined in (12.8p of a generic tensor 
A G cl”*’”']. Each tensor was generated using randn{n, ■ ■ ■ ,n) + i * randn{n, ■ ■ ■ ,n) in 
MATLAB. The results are reported in Table O in which N denotes the number of equiv¬ 
alence classes of isolated eigenpairs found by teig and T{m,n) denotes the bound of the 
number of equivalence classes of isolated eigenpairs (see Remark 12.4l( ai). 


(m, n) 

T{m, n) 

N 

time(s) 

(m, n) 

T{m, n) 

N 

time(s) 

(3,5) 

80 

80 

2.4 

(3,6) 

192 

192 

6.8 

(3,7) 

448 

448 

18.3 

(3,8) 

1024 

1024 

53.0 

(3,9) 

2304 

2304 

145.9 

(3,10) 

5120 

5120 

409.2 

(4,3) 

27 

27 

0.7 

(4,4) 

108 

108 

2.9 

(4,5) 

405 

405 

15.8 

(4,6) 

1458 

1458 

80.0 

(4,7) 

5103 

5103 

385.9 

(4,8) 

17496 

17496 

2115.5 

(5,3) 

48 

48 

1.2 

(5,4) 

256 

256 

8.8 

(5,5) 

1280 

1280 

73.8 

(5,6) 

6144 

6144 

606.5 

(5,7) 

28672 

28672 

5394.2 

(6,3) 

75 

75 

2.3 

(6,4) 

500 

500 

21.0 

(6,5) 

3125 

3125 

287.7 

(6,6) 

18750 

18750 

3721.3 

(7,3) 

108 

108 

3.6 

(7,4) 

864 

864 

51.5 

(7,5) 

6480 

6480 

981.3 


Table 2: Performance of teig on computing eigenpairs (|2.8I) of complex random tensors 


EXAMPLE 5.3 In this example we show the effectiveness and efficiency of eeig for 
finding all equivalence classes of isolated E-eigenpairs of a generic tensor A G Cl”*’”l. Each 
generic tensor was generated using randn{n, ■ ■ ■ , n)+i*randn{n, ■ ■ ■ , n) in MATLAB. The 
results are reported in Table [3l in which N denotes the number of equivalence classes of 
E-eigenpairs found by eeig and E{m,n) denotes the bound of the number of equivalence 
classes of isolated E-eigenpairs (see Remark l2.4l bB. 

According to m, [n, and [5], for a randomly generated tensor A G it has 

T{m,n) equivalence classes of eigenpairs (12.81) and E(m,n) equivalence classes of E- 
eigenpairs. Moreover, its eigenpairs and E-eigenpairs are isolated. From Tables [2] and 
[3l we observe that teig and eeig can find all equivalence classes of eigenpairs (12.81) and 
E-eigenpairs of such an tensor in the examples we tested. 

EXAMPLE 5.4 In this example we show the effectiveness and efficiency of teig and 
teneig for finding all equivalence classes of isolated R-eigenpairs of a tensor A, where 
A G C[”*’”1,R G are generic tensors. Each generic tensor was generated using 
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(m, n) 

E{m, n) 

N 

time(s) 

(m, n) 

E{m, n) 

N 

time(s) 

(3,5) 

31 

31 

1.4 

(3,6) 

63 

63 

3.1 

(3,7) 

127 

127 

7.5 

(3,8) 

255 

255 

20.3 

(3,9) 

511 

511 

48.5 

(3,10) 

1023 

1023 

133.9 

(4,3) 

13 

13 

0.4 

(4,4) 

40 

40 

1.7 

(4,5) 

121 

121 

5.4 

(4,6) 

364 

364 

26.9 

(4,7) 

1093 

1093 

119.5 

(4,8) 

3280 

3280 

555.8 

(5,3) 

21 

21 

0.7 

(5,4) 

85 

85 

4.2 

(5,5) 

341 

341 

22.3 

(5,6) 

1365 

1365 

166.5 

(5,7) 

5461 

5461 

1330.7 

(6,3) 

31 

31 

1.2 

(6,4) 

156 

156 

9.5 

(6,5) 

781 

781 

100.4 

(6,6) 

3906 

3906 

990.2 

(7,3) 

43 

43 

1.9 

(7,4) 

259 

259 

21.3 

(7,5) 

1555 

1555 

245.0 


Table 3: Performance of eeig on computing E-eigenpairs of complex random tensors 


randn{n, ■ ■ ■ ,n) + i * randn{n, ■ ■ ■ ,n) in MATLAB. The results are reported in Table 01 
in which N denotes the number of equivalence classes of eigenpairs found by teig or 
teneig, T{m,n) denotes the bound of the number of equivalence classes of isolated B- 
eigenpairs for m = m', and G{m,m',n) denotes the bound of the number of equivalence 
classes of isolated ;B-eigenpairs for m ^ m' (see Theorem 12.31) . 

From Tabled] we see that our teig and teneig find all equivalence classes of isolated 
0-eigenpairs of A for the generic tensors A and B we tested in a reasonable amount of 
time. 

5.2 Examples for Computing Real Eigenpairs 

In this subsection, numerical examples are provided to illustrate the effectiveness and 
efficiency of zeig or heig for computing real Z-eigenpairs or H-eigenpairs of a tensor 
A G By De fin ition 12.11 (A,x) is a Z-eigenpair if and only if (( ——x) is a 

Z-eigenpair, and (A, x) is an H-eigenpair if and only if (A, tx) is an H-eigenpair for any 
nonzero t G M. Only one representative from each equivalence class of eigenpairs will 
be listed in our examples. The notation A^^^ is used to denote I eigenvectors counting 
multiplicities are found for the eigenvalue A. In the following tables, the multiplicity of an 
eigenpair means the multiplicity of this eigenpair as a zero of the corresponding defining 
polynomial system. For the sake of conciseness, the polynomial system resulted from the 
tensor eigenvalue problem will be omitted. 

EXAMPLE 5.5 Consider the symmetric tensor A G whose corresponding polyno¬ 
mial form is the Motzkin polynomial: 

Ax^ = X^+ xfx2 + xfx2 — 3xfx2x|. 
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teig (m = m') 

teneig (m ^ m') 

{m, n) 

T(m, n) 

N 

time(s) 

(m, m', n) 

G{m, m', n) 

N 

time(s) 

(3,7) 

448 

448 

23.7 

(3,2,7) 

127 

127 

10.3 

(3,8) 

1024 

1024 

68.3 

(3,4,6) 

665 

665 

68.1 

(3,9) 

2304 

2304 

210.3 

(3,5,5) 

496 

496 

49.5 

(4,5) 

405 

405 

20.8 

(4,2,6) 

364 

364 

28.9 

(4,6) 

1458 

1458 

110.4 

(4,3,5) 

211 

211 

13.2 

(4,7) 

5103 

5103 

737.5 

(4,5,4) 

175 

175 

9.5 

(5,5) 

1280 

1280 

97.9 

(5,4,5) 

781 

781 

83.9 

(5,6) 

6144 

6144 

623.6 

(5,6,3) 

61 

61 

2.6 

(6,4) 

500 

500 

29.9 

(6,5,4) 

369 

369 

30.7 

(6,5) 

3125 

3125 

449.4 

(6,7,3) 

91 

91 

6.0 

(7,3) 

108 

108 

4.4 

(7,6,4) 

671 

671 

77.1 

(7,4) 

864 

864 

77.6 

(7,8,3) 

127 

127 

9.4 


Table 4: Performance of teig and teneig on computing generalized eigenpairs of complex 
random tensors 

In Example 5.9 of [5], it states that this tensor has 25 equivalence classes of Z-eigenpairs. 
Using zeig exactly 25 equivalence classes of Z-eigenpairs are found as shown in Table [5l 
which confirms the results of [5j. zeig takes about 0.9 seconds to carry out the entire 
computation. 


A 

0 ( 14 ) 

0.0156^^^ 

0.2500^^^ 

1 

Xl 

0.5774 

1 

0 

0.8253 

0.2623 

0.7071 

0 

X2 

±0.5774 

0 

1 

±0.2623 

±0.8253 

±0.7071 

0 

X3 

±0.5774 

0 

0 

±0.5000 

±0.5000 

0 

1 

multiplicity 

1 

5 

5 

1 

1 

1 

1 


Table 5: Z-eigenpairs of the tensor in Example 15.51 


All the H-eigenpairs found by Example 4.10 of [12] are also found by heig as shown 
in Table [H heig takes about 1.7 seconds to carry out the entire computation. 

So far the only available method that can find all real eigenvalues of a symmetric tensor 
is Algorithm 3.6 in m- In the next two examples, we report our experiments on zeig for 
computing all Z-eigenvalues, using examples taken from m- 

EXAMPLE 5.6 We use our zeig to compute the Z-eigenpairs of 12 symmetric tensors 
from |12] . The test problems and numerical results are given in the Appendix. From the 
numerical results we see that zeig finds all the Z-eigenvalues found by Algorithm 3.6 of 
m on this set of test problems. We now summarize the CPU time (in seconds) used by 
zeig in Table [T] Since the computer used in m is different from the computer used in 
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A 

0(i4) 

0 . 05551 ^^ 

l(15) 

Xl 

±1 

1 

0 

±0.4568 

1 

0 

1 

X2 

1 

0 

1 

1 

±0.4568 

0 

±1 

X3 

±1 

0 

0 

±0.6856 

±0.6856 

1 

1 

multiplicity 

1 

5 

5 

1 

1 

13 

1 


Table 6: H-eigenpairs of the tensor in Example 15.51 


this paper, the CPU time used by Algorithm 3.6 in m is not reported here, but we refer 

to [12]. 


Problem 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

time(s) 

0.3 

4.0 

0.3-0.4 

0.1 

0.6 

1.8 

15.7 

6.1 

0.3 

6.3 

27.3 

4.5 


Table 7: CPU time used by zeig for computing the Z-eigenvalues of 12 symmetric tensors 
from [T2| 


EXAMPLE 5.7 Consider the symmetric tensor A G (Example 4.16 in [T2|) with 

the polynomial form 

Ax‘^ = {xi - X 2 Y H-1- ( 2:1 - XnY + {X 2 - H-1- (X 2 - 

T ■ ■ ■ T (X)2_i Xn) • 

For different n, all the Z-eigenvalues found by Algorithm 3.6 in m are also found by 
zeig, which are given in Table [HI We remark that when n = 8,9,10, zeig can find all 
the Z-eigenvalues in a reasonable amount of time, but m reports that Algorithm 3.6 can 
only find the first three largest Z-eigenvalues. The CPU time used by zeig is reported in 
Table El Since different computers were used, we refer to m for the CPU time used by 
Algorithm 3.6 of [12] . For the sake of conciseness, the corresponding Z-eigenvectors are 
not displayed. 
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A 
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Appendix 

PROBLEM 1 Consider the symmetric tensor A € (Example 4.1 in [12], see also 

m) with the polynomial form 

Ax^ = x\ + 2 x 2 + 33^3 • 

Using zeig all the Z-eigenpairs found in [T2| are obtained (see Table[9|). zeig takes about 
0.3 seconds to carry out the entire computation. 


A 

0.5455^^^ 

0.6667^^^ 

0.7500^^^ 

1 

1.20 

2 

3 

Xl 

0.7385 

0.8165 

0.8660 

1 

0 

0 

0 

X2 

±0.5222 

±0.5774 

0 

0 

0.7746 

1 

0 

3^3 

±0.4264 

0 

±0.5000 

0 

±0.6325 

0 

1 

multiplicity 

1 

1 

1 

1 

1 

1 

1 


Table 9: Z-eigenpairs of the tensor in Problem [T] 


PROBLEM 2 For the diagonal tensor V G (Example 4.2 in [T2|) such that Vx^ = 
xf -|- 2^2 — Sxg — 4x|. Consider the symmetric tensor A G R^^d] guch that Ax^ = T>{Qx)^ 
where 

Q = {I — 2wiwJ){I — 2 w2W^){I — 2w3wJ) 

and wi,W 2 ,W 3 are randomly generated unit vectors. All the 30 Z-eigenpairs found in m 
are also found by using zeig. The 15 nonnegative Z-eigenvalues are listed below 

0.2518, 0.3261, 0.3466, 0.3887, 0.4805, 0.5402, 0.5550, 0.6057, 

0.8543, 0.9611, 1.0000, 1.2163, 2.0000, 3.0000, 4.0000. 

For conciseness, the corresponding Z-eigenvectors are not displayed here, zeig takes about 
4.0 seconds to do the entire computation. 

PROBLEM 3 Consider the symmetric tensor A G (Example 4.3 in [12], see also 

m) with the polynomial form 

Ax^ = 2xf + 3x2 + 5^3 + 4ax^X2X3, 

where a is a parameter. All Z-eigenvalues found in [12] are also found by zeig for different 
values of a, as shown in Table do] For conciseness, the corresponding Z-eigenvectors are 
not displayed here. The CPU time used by zeig for each a is also reported in the table. 

PROBLEM 4 Consider the symmetric tensor A G RI^’^1 (Example 4.4 in [12], see also 
|37j ) with the polynomial form 

Ax'^ = 3xf -I- X 2 + Qax\x\, 

where a is a parameter. All Z-eigenvalues found in m are also found by zeig for different 
values of a, which are listed in Table m The CPU time used by zeig for each a is also 
given in the table. For conciseness, the corresponding Z-eigenvectors are not displayed 
here. 
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a 

A 

time(s) 

0 

0.9677^), 1.20000, 1.42860, 1.87500, 2, 3, 5 

0.4 

0.25 

0.84640, 1.08810, 1.21500, 1.44120, 1.87500, 2, 3, 5 

0.4 

0.5 

0.72430, 1.20690, 1.25930, 1.47830, 1.87500 , 2, 3, 5 

0.4 

1 

0.47870, I. 6 I 33 O, 1.87500, 2 , 3, 5 

0.3 

3 

-0.51260, 1.87500, 2 , 2.21470, 3 , 5 

0.3 


Table 10: Z-eigenvalues of the tensor in Problem [3] 


a 

A 

time(s) 

-1 

-O. 6 OOOO, 1 , 3 

0.1 

0 

0.75000, 1 , 3 

0.1 

0.25 

0.97500, 1 , 3 

0.1 

0.5 

1, 3 

0.1 

2 

1, 3, 4.12500 

0.1 

3 

1, 3, 5 . 5714 O 

0.1 


Table 11: Z-eigenvalues of the tensor in Problem 0] 


PROBLEM 5 Consider the symmetric tensor A G (Example 4.5 in [12], see also 

[23] or [34]) such that 

Ann = 0.2883, Ani2 = -0.0031, Ania = 0.1973, An22 = -0.2485, 

^1123 = —0.2939,^1133 = 0.3847, Ai222 = 0.29 72, A 1223 = 0.1862, 

^1233 = 0.091 9, A 1333 = —0.3619,^2222 = 0.1241,^2223 = —0.34 20 , 

A 2233 = 0.2127,^2333 = 0.27 27,^3333 = —0.30 54. 

All the Z-eigenpairs found in m are also found by zeig, as given in Table [12] zeig takes 
about 0.6 seconds to do the entire computation. 


A 

-1.0954 

-0.5629 

-0.0451 

0.1735 

0.2433 

0.2628 

0.2682 

0.3633 

0.5105 

0.8169 

0.8893 

Xi 

-0.5915 

-0.1762 

0.7797 

0.3357 

-0.9895 

-0.1318 

0.6099 

0.2676 

-0.3598 

-0.8412 

-0.6672 

X2 

0.7467 

0.1796 

0.6135 

0.9073 

-0.0947 

0.4425 

0.4362 

0.6447 

0.7780 

0.2635 

-0.2471 

^3 

0.3043 

-0.9678 

0.1250 

0.2531 

0.1088 

0.8870 

0.6616 

0.7160 

-0.5150 

-0.4722 

0.7027 


Table 12: Z-eigenpairs of the tensor in Problem [5] 


PROBLEM 6 Consider the symmetric tensor A G (Example 4.6 in [12], see also 

[39]) such that Am = z for i = 1 ,..., 6 and = 10 for z = 1,...,5 and zero 

otherwise. All the Z-eigenvalues found in [12] are also found by zeig. The 19 nonnegative 
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Z-eigenvalues are listed below: 

3.9992 4.0225 4.2464 4.3358 5.1402 5.4817 5.5218 5.5668 

5.5674 6.0000 7.2165 8.1889 8.5979 8.6596 8.7347 10.9711 

15.4298 15.4552 16.2345 

For conciseness, the corresponding Z-eigenvectors are not displayed here, zeig takes about 
1.8 seconds to carry out the entire computation. 

PROBLEM 7 Consider the symmetric tensor A G (Example 4.7 in [12], see also 

|28j ) with the polynomial form 

-Ax^ = (Xl - X2)^ + (Xl - + (Xl - X4)^ + (xi - xs)^ + (xi - Xg)^ 

+ (X2 - Xa)"^ + (X2 - X4)^ + (X2 - Xs)^ + (X2 - Xg)"^ 

+(X3 - X4)"^ + (xa - Xs)"^ + (xa - xg)^ 

+ (X4 - Xsf + (X4 - Xg)'^ + (X5 - Xg)^. 

All the 5 Z-eigenvalues found in [T2| are also found by zeig, which are given in Table [T3j 
As pointed out in [T2|, every permutation of a Z-eigenvector is also a Z-eigenvector. Only 


A 

x^ 

multiplicity 

-7.2000^^^^ 

(0.1826,0.1826,0.1826,0.1826, 0.1826, -0.9129) 

1 

- 6 . 0000 ^ 1 ^^ 

(0.7071,0,0,0,0,-0.7071) 

1 

-4.5000^ 

(0.5774,0.5774, -0.2887, -0.2887, -0.2887, -0.2887) 

- 

-4.0000^1*^^ 

(0.4082, 0.4082, 0.4082, -0.4082, -0.4082, -0.4082) 

1 


(0.4082, 0.4082, 0.4082,0.4082,0.4082,0.4082) 

- 


Table 13; Z-eigenpairs of the tensor in Problem [7] 


the Z-eigenvector with xi > X 2 > • • • > xg corresponding to one Z-eigenvalue is listed. 
We remark that the Z-eigenpairs corresponding to Z-eigenvalues 0 and —4.5 are in a posi¬ 
tive dimensional solution component of the corresponding polynomial system. Therefore, 
there are infinitely many Z-eigenvectors associated with 0 and —4.5. zeig finds 484 Z- 
eigenvectors associated with 0 and 180 Z-eigenvectors associated with —4.5. Only one of 
these Z-eigenvectors for each case is listed in Table (Hi zeig takes about 15.7 seconds to 
do the entire computation. 

PROBLEM 8 Consider the symmetric tensor A G Rt^’^1 (Example 4.8 in [T2|, see also 
|46j i with the polynomial form 

Ax^ = (xi -I- X 2 + Xa + X 4 )'‘ -I- (x 2 + Xa + X 4 -I- xg)"^. 

All the 3 Z-eigenvalues found in [T 2 | are also found by zeig, which are shown in Table [TTl 
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A 

x^ 

multiplicity 

0^*1 

(0.3870, -0.1537,0.4532, -0.6866,0.3870) 

- 

0.5000 

(0.7071,0,0,0,-0.7071) 

1 

24.5000 

(0.2673,0.5345,0.5345,0.5345,0.2673) 

1 


Table 14: Z-eigenpairs of the tensor in Problem [ 8 ] 


We remark that the Z-eigenpairs corresponding to Z-eigenvalue 0 are in a positive 
dimensional solution component of the corresponding polynomial system. Thus, there are 
infinitely many Z-eigenvectors associated with Z-eigenvalue 0. zeig finds 234 of them. 
Only one of them is listed in Table [TH zeig uses about 6.1 seconds to do the entire 
computation. 

PROBLEM 9 Consider the symmetric tensor A G (Example 4.9 in [T^], see also 

0 ) with the polynomial form 

Ax — -j- 3x^X2 ~h 3xix^. 

We remark that the Z-eigenpairs corresponding to Z-eigenvalue 2 are in a positive di¬ 
mensional solution component of the corresponding polynomial system. Thus, there are 
infinitely many Z-eigenvectors associated with Z-eigenvalue 0. zeig finds 7 of them. Only 
one of them is listed in Table [T^ z:eig uses about 0.3 seconds to do the entire computation. 


A 

x^ 

multiplicity 

2 W 

( 1 , 0 , 0 ) 

- 


Table 15; Z-eigenpairs of the tensor in Problem [9] 


PROBLEM 10 Consider the tensor A € (Example 4.12 in [12], see also [M]) such 
that 

Ai,..,j4 = sin(zi +i2 + i3 + k)- 

When n = 5, all the 5 Z-eigenvalues found in [12] are also found by zeig, which are given 
in Table dU 

We remark that the Z-eigenpairs corresponding to Z-eigenvalue 0 are in a positive 
dimensional solution component of the corresponding polynomial system. Thus, there are 
infinitely many Z-eigenvectors associated with 0. zeig finds 234 of them. Only one of them 
is listed in Table fTHl zeig takes about 6.3 seconds to carry out the entire computation. 

PROBLEM 11 Consider the tensor A G (Example 4.13 in [T^]) such that 

= tan(zi) -I-|-tan(i 4 ). 
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A 


multiplicity 

-8.8463 

(0.5809,0.3563, -0.1959, -0.5680, -0.4179) 

1 

-3.9204 

(-0.1785,0.4847,0.7023, 0.2742, -0.4060) 

1 

QP) 

(-0.5213,0.3748, -0.6608,0.1824, -0.3433) 

- 

4.6408 

(0.5055, -0.1228, -0.6382, -0.5669,0.0256) 

1 

7.2595 

(0.2686,0.6150,0.3959, -0.1872, -0.5982) 

1 


Table 16: Z-eigenpairs of the tensor in Problem [TUI 


A 


multiplicity 

-133.2871 

(0.1936, 0.5222, 0.3429,0.2287,0.6272,0.3559) 

1 

OP) 

(-0.5840, -0.3454,0.1784, 0.6773,0.1892, -0.1156) 

- 

45.5045 

(0.6281,0.0717,0.3754,0.5687, -0.1060, 0.3533) 

1 


Table 17: Z-eigenpairs of the tensor in Problem [TT] 


When n = 6, all the 3 Z-eigenvalues found in [12] are also found by zeig, which are given 
in Table fTTl 

We remark that the Z-eigenpairs corresponding to Z-eigenvalue 0 are in a positive 
dimensional solution component of the corresponding polynomial system. Thus, there 
are infinitely many Z-eigenvectors associated with 0. zeig hnds 724 of them. Only one 
of them is listed in Table [T71 It takes zeig about 27.3 seconds to carry out the entire 
computation. 

PROBLEM 12 Consider the tensor A € (Example 4.14 in [T2|) such that 

Ai,...,i5 = In(ii) H-h ln(i5). 

For n = 4, all the 3 Z-eigenvalues found in [T2| are also found by zeig, which are shown 
in Table mi 


A 

x^ 

multiplicity 

0^*1 

(-0.4304,0.8139,0.0069, -0.3903) 

- 

0.7074 

(-0.9054, -0.3082,0.0411,0.2890) 

1 

132.3070 

(0.4040, 0.4844, 0.5319,0.5657) 

1 


Table 18: Z-eigenpairs of the tensor in Problem [12] 


We remark that the Z-eigenpairs corresponding to Z-eigenvalue 0 are in a positive 
dimensional solution component of the corresponding polynomial system. Thus, there are 
infinitely many Z-eigenvectors associated with 0. zeig finds 166 of them. Only one of 
them is listed in Table mi The entire computation takes zeig about 4.5 seconds. 


35 






















